#include "calc.h"

// Compilation Routines
Calculator::Calculator()
{
	newLine=FALSE;
	strDepth=3;
	defUserf=NULL;
	defVar=NULL;
	defUserfNames=NULL;
	defVarNames=NULL;
	first=NULL;
	last=NULL;
	input=NULL;
	count=0;
}
Calculator::~Calculator()
{
#ifdef DEBUG_MATH
	int B=count*sizeof(Expression);
	int KB=B/1024;
	B%=1024;
	int MB=KB/1024;
	KB%=1024;
	_tprintf(_T("Expressions To Clear: %d (%dMB %dKB %dB)\n"),count,MB,KB,B);
#endif
	clear();
}
void Calculator::clear() // delete All Expressions
{
	if(!first)
		return;
	Expression *expr=first;
	while(expr->next)
	{
		expr=expr->next;
		delete expr->prev;
	}
	delete expr;
	count=0;
	first=NULL;
	last=NULL;
}
void Calculator::del(Expression * expr)
{
	if(!expr)
		return;
	if(last==expr)
		last=expr->prev;
	if(first==expr)
		first=expr->next;
	count--;
	if(expr->prev)
		expr->prev->next=expr->next;
	if(expr->next)
		expr->next->prev=expr->prev;
	delete expr;
}
// Create An Empty Expression (0)
Expression * Calculator::create()
{
	Expression * expr=new Expression();
	if(!expr)
		return NULL;
	expr->prev=last;
	if(last!=NULL)
		last->next=expr;
	if(first==NULL)
		first=expr;
	last=expr;
	count++;
#ifdef DEBUG_MATH
	expr->id=count;
#endif
	return expr;
}
// expr = (type!=0) factor.A + cnst, (type=0) factor.A
// A = B.mul + add
// B = coef.func(C)^pow
// C = (type<0) child, (type>0) var, (type=0) cnst
Expression * Calculator::create(
					int type,
					int func,
					double coef,
					double pow,
					double cnst,
					double factor,
					Expression *add,
					Expression *mul,
					Expression *child)
{
	Expression * expr=create();
	if(expr)
	{
		expr->add=add;
		expr->mul=mul;
		expr->child=child;

		expr->coef=(float)coef;
		expr->factor=(float)factor;
		expr->pow=(float)pow;
		expr->cnst=(float)cnst;

		expr->func=func;
		expr->type=type;
	}
	return expr;
}
// Create And Copy
Expression * Calculator::create(const Expression * model)
{
	if(!model)
		return NULL;
	Expression * expr=create();
	if(expr)
	{
		expr->add=model->add;
		expr->mul=model->mul;
		expr->child=model->child;

		expr->coef=model->coef;
		expr->factor=model->factor;

		expr->cnst=model->cnst;

		expr->func=model->func;
		expr->pow=model->pow;
		expr->type=model->type;

		expr->d=model->d;
	}
	return expr;
}
// Shortcut To Create Constant Cnst
Expression * Calculator::createc(float cnst)
{
	Expression * expr=create();
	if(expr)
	{
		expr->cnst=cnst;
	}
	return expr;
}
// Shortcut To Create Parent Expression (0)
Expression * Calculator::createp(Expression *child, double pow, BOOL ifNeeded, BOOL trueChild)
{
	if(!child)
		return NULL;
	if(ifNeeded&&!needParent(child))
		return child;
	if(trueChild)
		child=getTrueChild(child);
	Expression * expr=create();
	if(expr)
	{
		expr->type=-1;
		expr->child=child;
		expr->pow=(float)pow;
	}
	return expr;
}
// Shortcut To Create Variable 0
Expression * Calculator::createv(int varID, double coef, double pow)
{
	if(varID<=0)
		return NULL;
	Expression * expr=create();
	if(expr)
	{
		expr->type=varID;
		expr->coef=(float)coef;
		expr->pow=(float)pow;
	}
	return expr;
}
// Shortcut To Create A Function f(0)
Expression * Calculator::createf(int funcID, int type, Expression *child, int diffs)
{
	if(!funcID)
		return NULL;
	Expression * expr=create();
	if(expr)
	{
		expr->type=type;
		expr->func=funcID;
		expr->child=child;
		expr->d=diffs;
	}
	return expr;
}
double Calculator::evalFunc(const Expression *expr, double x)
{
	if(expr->func<0&&userf)
	{
		dfdx f=userf[-1-expr->func];
		if(f)
			return f(x,expr->d);
	}
	else switch(expr->func)
	{
	case FUNC_EXP:
		return exp(x);
	case FUNC_LOG:
		return log(x);
	case FUNC_SIN:
		return sin(x);
	case FUNC_COS:
		return cos(x);
	case FUNC_TAN:
		return tan(x);
	case FUNC_ASIN:
		return asin(x);
	case FUNC_ACOS:
		return acos(x);
	case FUNC_ATAN:
		return atan(x);
	}
	return x;
}
void Calculator::addFuncName(DynamicString& str, const Expression *expr)
{
	if(expr->func<0)
	{
		if(userfNames)
		{
			str+=userfNames[-expr->func-1];
		}
		else
		{
			str+='f';
			str+=format.toString(int(-expr->func));
		}
		if(expr->d)
		{
			str+='\'';
			str+=format.toString((int)expr->d);
		}
	}
	else switch(expr->func)
	{
	case FUNC_EXP:
		str+=_T("exp");
		break;
	case FUNC_LOG:
		str+=_T("log");
		break;
	case FUNC_SIN:
		str+=_T("sin");
		break;
	case FUNC_COS:
		str+=_T("cos");
		break;
	case FUNC_TAN:
		str+=_T("tan");
		break;
	case FUNC_ASIN:
		str+=_T("asin");
		break;
	case FUNC_ACOS:
		str+=_T("acos");
		break;
	case FUNC_ATAN:
		str+=_T("atan");
		break;
	}
}
// Create A True Func Expression, MUL = 0, ADD = 0, CNST = 0
Expression * Calculator::diffFunc(const Expression *expr)
{
	Expression *dfun;
	switch(expr->func)
	{
	case FUNC_NONE: // 1
		return createc(1);
	case FUNC_EXP: // exp(x)
		return createf(FUNC_EXP,expr->type,expr->child);
	case FUNC_LOG: // 1/x
		dfun=create();
		dfun->type=expr->type;
		dfun->child=expr->child;
		dfun->pow=-1;
		return dfun;
	case FUNC_SIN: // cos(x)
		return createf(FUNC_COS,expr->type,expr->child);
	case FUNC_COS: // -sin(x)
		dfun=createf(FUNC_SIN,expr->type,expr->child);
		dfun->coef=-1;
		return dfun;
	case FUNC_TAN: // 1/cos(x)^2
		dfun=createf(FUNC_COS,expr->type,expr->child);
		dfun->pow=-2;
		return dfun;
	case FUNC_ASIN: // 1/sqrt(1-x^2)
		dfun=create();
		dfun->type=expr->type;
		dfun->child=expr->child;
		dfun->pow=2;
		dfun->coef=-1;
		dfun->cnst=1;
		dfun=createp(dfun,-0.5);
		return dfun;
	case FUNC_ACOS:  // -1/sqrt(1-x^2)
		dfun=create();
		dfun->type=expr->type;
		dfun->child=expr->child;
		dfun->pow=2;
		dfun->coef=-1;
		dfun->cnst=1;
		dfun=createp(dfun,-0.5);
		dfun->coef=-1;
		return dfun;
	case FUNC_ATAN:  // 1/(1+x^2)
		dfun=create();
		dfun->type=expr->type;
		dfun->child=expr->child;
		dfun->pow=2;
		dfun->cnst=1;
		dfun=createp(dfun,-1);
		return dfun;
	}
	//Not Found, So Must Be User Defined
	return createf(expr->func,expr->type,expr->child,expr->d+1);
}
double Calculator::eval(Expression *expr, double *var, dfdx *userf) // Calculate Expression
{
	if(expr==NULL)
		return 0;
	input=expr;
	if(var==NULL)
		var=defVar;
	if(userf==NULL)
		userf=defUserf;
	this->var=var;
	this->userf=userf;
	return getValue(expr);
}
double Calculator::getValue(const Expression *expr)
{
	if(!expr->factor)
		return 0;
	double value=0;
	if(expr->coef)
	{
		if(expr->type)
		{
			if(expr->type<0&&expr->child!=NULL)
			{
				value=getValue(expr->child);
			}
			else if(expr->type>0&&var)
			{
				value=var[expr->type-1];
				value=evalFunc(expr,value);
			}
		}
		else
		{
			value=expr->cnst;
		}
		if(value)
		{
			if(expr->pow!=1&&value!=1)
				value=pow(value,(double)expr->pow);
			value*=expr->coef;
			if(expr->mul!=NULL)
			{
				value*=getValue(expr->mul);
			}
		}
	}
	if(expr->add!=NULL)
	{
		value+=getValue(expr->add);
	}
	if(expr->type)
		value+=expr->cnst;
	return expr->factor*value;
}
STR  Calculator::toString(const Expression *expr, STR *varNames,  STR *userfNames, int depth)
{
	if(depth==0)
		depth=this->strDepth;
	if(expr==NULL)
		return NULL;
	if(varNames==NULL)
		varNames=defVarNames;
	if(userfNames==NULL)
		userfNames=defUserfNames;
	this->userfNames=userfNames;
	this->varNames=varNames;
	dstr.cut();
	addStr(dstr,depth,expr);
	str.create(dstr.length());
	dstr.toString(str);
	return str;
}
void Calculator::addStr(DynamicString& str, int depth, const Expression *expr)
{
	BOOL zero=TRUE, parent;
	if(!expr->factor)
	{
		str+='0';
		return;
	}
	parent=(expr->factor!=1)||needParent(expr);
	if(parent)
	{
		if(expr->factor==-1)
		{
			str+='-';
		}
		else if(expr->factor!=1)
		{
			if(expr->factor==floor(expr->factor))
				str+=format.toString((int)expr->factor);
			else
				str+=format.toString(expr->factor);
		}
		str+='(';
	}
	if(expr->coef)
	{
		zero=FALSE;
		if(expr->coef!=1)
		{
			if(expr->coef==-1)
				str+='-';
			else
			{
				if(expr->coef==floor(expr->coef))
					str+=format.toString((int)expr->coef);
				else
					str+=format.toString(expr->coef);
			}
		}
		if(expr->func)
		{
			addFuncName(str,expr);
			str+='(';
		}
		if(expr->type)
		{
			if(expr->type<0&&expr->child!=NULL)
			{
				if(depth)
				{
					if(!expr->func)
						str+='(';
					addStr(str,depth-1,expr->child);
					if(!expr->func)
						str+=')';
				}
				else
					str+='A';
			}
			else
			{
				if(varNames)
				{
					str+=varNames[expr->type-1];
				}
				else
				{
					str+='x';
					str+=format.toString((int)expr->type);
				}
			}
		}
		else
		{
			if(!expr->func)
				str+='(';
			if(expr->cnst==floor(expr->cnst))
				str+=format.toString((int)expr->cnst);
			else
				str+=format.toString(expr->cnst);
			if(!expr->func)
				str+=')';
		}
		if(expr->func)
		{
			str+=')';
		}
		if(expr->pow!=1)
		{
			str+='^';
			if(expr->pow==floor(expr->pow))
				str+=format.toString((int)expr->pow);
			else
				str+=format.toString(expr->pow);
		}
		if(expr->mul)
		{
			str+='(';
			addStr(str,depth,expr->mul);
			str+=')';
		}
	}
	if(expr->add)
	{
		if(!zero&&expr->add->factor>0)
			str+='+';
		zero=FALSE;
		if(newLine)
			str+='\n';
		addStr(str,depth,expr->add);
	}
	if(expr->type&&expr->cnst)
	{
		if(!zero&&expr->cnst>0)
			str+='+';
		if(expr->cnst==floor(expr->cnst))
			str+=format.toString((int)expr->cnst);
		else
			str+=format.toString(expr->cnst);
	}
	if(parent)
	{
		str+=')';
	}
}
Expression * Calculator::diff(Expression *expr, int varID) // Create Derivate Expression
{
	if(expr==NULL&&varID<=0)
		return NULL;
	input=expr;
	this->varID=varID;
	return getDiff(expr);
}
int Calculator::clearUnused(Expression *expr)
{
	if(!expr)
		return 0;
	Expression *tmp=first;
	int count=0;
	while(tmp)
	{
		tmp->connect=NULL;
		tmp=tmp->next;
	}
	connectTree(expr);
	tmp=first;
	while(tmp->next)
	{
		tmp=tmp->next;
		if(!tmp->prev->connect&&tmp->prev!=expr)
		{
			del(tmp->prev);
			count++;
		}
	}
	if(!tmp->connect&&tmp!=expr)
	{
		del(tmp);
		count++;
	}
#ifdef DEBUG_MATH
	_tprintf(_T("Unused expressions removed: %d\n"),count);
#endif
	return count;
}
void Calculator::remExtraParent(Expression *expr)
{
	if(!expr)
		return;
	if(expr->add)
	{
		remExtraParent(expr->add);
	}
	if(expr->mul)
	{
		remExtraParent(expr->add);
	}
	if(expr->child)
	{
		expr->child=getTrueChild(expr->child);
		remExtraParent(expr->child);
	}
}
void Calculator::connectTree(Expression *expr)
{
	if(!expr)
		return;
	if(expr->add&&!expr->add->connect)
	{
		expr->add->connect=expr;
		connectTree(expr->add);
	}
	if(expr->mul&&!expr->mul->connect)
	{
		expr->mul->connect=expr;
		connectTree(expr->mul);
	}
	if(expr->child&&!expr->child->connect)
	{
		expr->child->connect=expr;
		connectTree(expr->child);
	}
}
Expression * Calculator::getDiff(const Expression *expr)
{	// expr' = factor. A'
	// A' = B'.mul + mul'.B + add'
	// B' = (coef.pow)func(C)^(pow-1).func'(C).C'
	// C' = (type<0) child', (type>0) 1, (type=0) 0
	Expression *B,*dA=NULL,*dadd=NULL,*dmul=NULL,*dB=NULL,*dchild=NULL,*dfunc=NULL;
	if(!expr||!expr->factor)
		return NULL;
	if(expr->add)
	{
		dA=dadd=createp(getDiff(expr->add),1,TRUE);
	}
	if(expr->coef)
	{
		if(expr->mul)
		{
			dmul=createp(getDiff(expr->mul),1,TRUE);
			if(dmul)
			{
				// B = coef.func(0)^pow
				B=create(expr);
				if(B->type)
					B->cnst=0;
				B->factor=1;
				B->add=NULL;
				B->mul=NULL;
				dmul->mul=B;
				dmul->add=dA;
				dA=dmul;
			}
		}
		// dB=coef.pow.func(0)^(pow-1).func'(0).0'
		if(expr->pow&&expr->type) // If Not Const
		{
			BOOL arg=TRUE;
			// 0'
			if(expr->type<0)
			{
				if(expr->child)
					dchild=createp(getDiff(expr->child),1,TRUE);
				if(!dchild)
					arg=FALSE;
			}
			else
			{
				if(expr->type!=varID)
					arg=FALSE;
			}
			if(arg) // If 0' Is Not Zero
			{
				if(expr->func)
				{
					// func'(0)
					dfunc=diffFunc(expr);
				}
				if(expr->pow!=1)
				{
					// (coef.pow).func(0)^(pow-1)
					dB=create();
					dB->type=expr->type;
					dB->func=expr->func;
					dB->pow=expr->pow-1;
					dB->child=expr->child;
				}
				if(dfunc)
				{
					dfunc->mul=dB;
					dB=dfunc;
				}
				if(dchild)
				{
					dchild->mul=dB;
					dB=dchild;
				}
				dB=createp(dB,1,TRUE);
				if(!dB)
					dB=createc(1);
				dB->coef*=expr->coef*expr->pow;
				dB->mul=expr->mul;
				dB->add=dA;
				dA=dB;
			}
		}
	}
	if(dA)
		dA->factor*=expr->factor;
	return dA;
}

int main(int argn, char ** argv)
{
	INITRAND;

	int i=0,j=0,k=0;
	FILE *f=NULL;
	_TCHAR fileName[128];
	PrecisionTimer t;

	const int N=12;
	int count,size;
	double maxDiffT=0,maxEvalT=0;
	STR str=NULL;
	STR varNames[2], userfNames[2];
	Calculator calc;
	const int _X=0,_E=1,_SIN=0;
	dfdx userf[1];
	double vars[2],nmin=-5,nmax=5,emin=0,emax=1,n=-1.5;
	double tmp;
	Expression *expr=NULL,*expr2=NULL;

	varNames[_X]=_T("x");
	varNames[_E]=_T("e");
	userfNames[_SIN]=_T("f");
	vars[_X]=PI/6;
	userf[_SIN]=df;
	calc.defVar=vars;
	calc.defUserf=userf;
	calc.defVarNames=varNames;
	calc.defUserfNames=userfNames;
	size=sizeof(Expression);
	_tprintf(_T("diff depth %d, expr size=%d\n"),N-1,size);
	vars[_E]=e;
	_stprintf(fileName,_T("diff%d.txt"),N);
	f=_tfopen(fileName,_T("w"));
	if(f)
	{
		_tprintf(_T("Writing values in file '%s'...\n"),fileName);
		_ftprintf(f,_T("Diff values for n=%.3f, e=%.8f :\n"),n,e);
	}
	expr=calc.create( 2,-1,1,n);
	int factN=1;
	for(i=1;i<=N;i++)
	{
		factN*=i;
		_tprintf(_T("Searching %2d derivate for n=%.3f... "),i,n);
		
		t.start();
		expr=calc.diff(expr,2);
		t.stop();
		count=calc.getExprCount();
		//str=calc.toString(expr);
		_tprintf(_T("DONE (%d exprs, time=%.3f ms)\n"),count,1e3*t.time);
		if(maxDiffT>0&&t.time>maxDiffT)
		{
			_tprintf(_T("Max Diff Time Exceeded! (exprs=%d)\n"),count);
			return 0;
		}
		_tprintf(_T("Evaluating %2d derivate for e=%.8f... "),i,e);
		t.stop();
		tmp=calc.eval(expr);
		calc.clearUnused(expr);
		t.stop();
		_tprintf(_T("DONE: value=%.16f, time=%.3f ms\n"),tmp,1e3*t.time);
		if(maxEvalT>0&&t.time>maxEvalT)
		{
			_tprintf(_T("Max Eval Time Exceeded!\n"),count);
			return 0;
		}
		if(f)
		{
			_ftprintf(f,_T("\ndiff%02d = %+.16e\n"),i,tmp);
			_ftprintf(f,_T("dc%02d = %+.16e/%d = %+.16e\n"),i,tmp,factN,tmp/factN);
			_ftprintf(f,_T("ic%02d = %+.16e/%d = %+.16e\n\n"),i+1,tmp,factN*(i+1),tmp/(factN*(i+1)));
		}
		//_tprintf(_T("Press ENTER to continue");
		//getchar();
	}
	if(f)
	{
		fclose(f);
		f=NULL;
	}

	return 0;
}
